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1. Introduction 

Since Wigner's seminal work of the 1950's, random matrix theory (RMT) has established 
itself as a cornerstone of modern theoretical physics, with innumerable applications 
(see, for example, [1] and references therein). One central problem of RMT is the 
determination of the mean spectral density of an ensemble of random matrices. In some 
cases, well understood universal laws governing the mean spectral density have been 
known for some time (see, for instance, [21 EJ H]). However, for ensembles of sparse 
matrices, i.e. matrices with many entries being zero, the picture is rather different. 
Studied first by Rodgers and Bray [5], the spectral density of real symmetric sparse 
random matrices has been extensively researched [6[ [7J [8j [9[ [TOj [TTJ [12] , although exact 
results have been only obtained relatively recently [131 HH HHJ EE fTTj . 
One area of research in which sparse random matrices feature heavily is in the modelling 
of real-world complex networks, with applications in fields as diverse as bioinformatics 
and finance. In the investigation of spectral density, the most commonly studied 
random graph ensembles are those in which the degrees of neighbouring vertices become 
independent in the limit iV — ► oo. This class includes the classical Erdos-Reyni (or 
Poissonian) random graphs (5j EJ [7J [El EEU [El EE] , and graphs with a specified degree 
distribution [TH [17]. Unfortunately, these simple ensembles may not provide realistic 
models for real-world complex networks which can include features such as correlations 
between non-neighbouring vertices, or groups of vertices organised into highly connected 
communities [15] . Very few results have been obtained for random graph ensembles with 
more complex topologies, though there exists an approximation scheme for graphs with 
degree-degree correlations [8] and numerical investigations of some other ensembles [19J. 
In this paper, we extend the analysis of the spectral density of random graphs to complex 
networks with the features described above. In order to achieve this we consider graphs 
with hierarchically constrained topologies. Introduced in [2D1 El], such ensembles can 
be tuned to more closely reflect the statistics of a real graph, whilst remaining in a form 
for which large N computations are tractable. We also propose a simple generalisation 
of this ensemble to one featuring a community structure. 

Applying the replica method, expressions are found for spectral densities of these 
ensembles in terms of the statistical properties of the graphs. For the sake of clarity we 
only show here the calculations for the spectral density of connectivity matrices, however 
the techniques can just as well be applied to more general sparse random matrices. This 
is possible even in the case of non-Hermitian matrices, following a scheme similar to 
that set out in [22J. 

The ensemble definitions are given in Section 2. In the third section, we prepare for the 
calculation of the spectral density by first computing relevant statistical properties of 
the graph ensembles, gaining some insight to the problem on the way. Section 4 contains 
the replica calculation itself, both for the hierarchically constrained ensemble and for 
graphs with a community structure. Several applications are discussed in Section 5, 
where analytical findings of the resulting spectral densities are compared with numerical 
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diagonalisations. The final section contains a summary and discussion. 
2. Ensemble definitions 

Let G = (V, E) be a graph on a set V = {1, N} of vertices and E C V x V of edges. 
The latter set is usually represented by the connectivity matrix C, whose entries are 
Cij = 1 if (i, j) G -E, and Gy = otherwise Vi, j € V. We define the £-th generalised 
degree of a vertex i, denoted as kf (C), to be the total number of walks of length i 
starting at vertex i. In terms of the connectivity matrix, this is given by the recursive 
definition 

JV 

*f(C) = l, kf\C) = J2^kf~ 1 \C), £=1,...,L, (1) 

for some integer L. We denote ki(C) = (kf (C), ■ ■ ■ , k\ L \C)) . Note that the first 
component of fcj(C) is simply the degree of vertex i, so we will usually refer to it 
as ki(C), rather than Kp(C). The generalised degree of a particular vertex contains 
information about the generalised degrees of its neighbours. For instance, if vertex i has 
degree k and neighbours {ji, jk}, then 

fcf ) (C) = E^r i) (C) for£=l,...,L. (2) 
t=i 

In what follows we study graph ensembles in which the generalised degrees are 
constrained. To achieve this, we follow the scheme set out in [201 12T] and define the 
weight of the graph ensemble as 

N 

W N (C) = n [^Q{K h,)^ A + (l - ^Q{k u k 3 )) 8 Cij , ] J] 8 kiio) , ki , (3) 

i<j i=l 

where the {ki)f =1 are taken to be arbitrary, c is the average c = (1/N) ^2i=x hi an( i 
Q(k, k') is a symmetric, non-negative function. Note from ([3]) that for each vertex i, its 
generalised degree fcj(C) is constrained to be precisely fej. For this reason we also refer 
to the {ki}^ as generalised degrees. 

As mentioned in the introduction, this ensemble is designed to be amenable to analysis 
in the large iV limit, whilst allowing the topology of the graphs to be tuned by choosing 
the generalised degrees {fcj} and the function Q(k,k'). This was demonstrated in [21] 
with the computation of the entropy of the ensemble, and in [23] the Ising model is 
analysed on such graphs in the simpler case of L = 1. 

We will also consider an ensemble whose graphs feature community structures. In such 
graphs, vertices are organised into densely intra-connected clusters, also called modules 
or communities, with a sparse distribution of inter-community edges 0- To incorporate 
this structure without sacrificing the solvability of the model, we propose a generalisation 

| For a nice and complete review on community structures and its importance in complex networks 
see [18]. 
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of the previously introduced ensemble %. 

Consider a graph composed of iV communities, each of size M (giving a total of NM 
vertices) . We decompose the connectivity matrix C of this graph into three sets of smaller 
matrices: a single N x N connectivity matrix C with entries qj = 1 if communities i 
and j are connected, and zero otherwise; a collection ofMxM matrices Bij encoding 
the connections between vertices in communities % and j\ and a collection of M x M 
matrices Aj specifying the internal connections of community i. For our ensemble, we 
take the B^ and A; to be drawn randomly and independently according to weights 
n{B) and u(A), respectively, whilst C is taken from the constrained generalised degree 
ensemble. All together, the weight for this graph ensemble can be written as follows 

N 

Wfi%(C) = J] [jjQ^ k^.MBij) + (l - kjj) S^jo] n v(Ai)5 k .(c),ki (4) 

i<j i=l 

For ease of calculation, we take ji to satisfy f-i(B) = n(B T ). 

Exploiting the bridge to statistical mechanics introduced by Edwards and Jones |27J, 
we will compute the mean spectral density of the connectivity matrices of graphs from 
both ensembles in the limit A" — > oo by the replica method. It will be instructive to 
first analyse the local statistics of the graphs in the same limit, as this will help us to 
simplify the subsequent analysis. 



3. Asymptotic graph statistics 

For the ensembles of graphs under study, the relevant statistical properties are captured 
by the following joint distribution 

P({Qt}!Li,k) =^ (^2 <J **(°).* 6 kn(o), qi ---^ k (C), qk c ijl ---c ijk \ , (5) 

\ »=i h<-<ik I c 

where here, and hereafter, we use (• • -)c to denote the ensemble average. Note that 
(jHJ) is the probability of finding a vertex whose neighbour has generalised degree k and 
is connected to k vertices with generalised degrees {qt}t=i- Other relevant probability 
distributions can be obtained by marginalising the expression (jSJ). In particular, by 
summing with respect to {q t }f =1 and denoting the resulting distribution as P(k) we 
obtain 

p ( fe )=i™(^EW)) . (6) 

\ i=l I c 

which is the probability of finding a vertex connected to a vertex with generalised 
degree k (see, for instance, [21]). Note that since in our ensemble the generalised 
degrees are constrained to the arbitrary values {fcij^Lx, we have that P(k) = p(k)k/c 
with p(k) = hmjv->oo(l/-^0 Si=i being the generalised degree distribution. Other 



§ Very recently, the spectral density of matrices with a particular type of modular form has been 
studied in [H[25l|26]. 
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probability distributions, however, may not have such straightforward forms. 
An expression for the distribution (jSJ) can be found through a saddle-point computation, 
the details of which will prove to be of great use in the later replica analysis of the 
spectral density, as our final equations for the spectral density will be written in terms 
of certain marginals of P ({qt\t=ii 



3.1. General calculation 



To find an expression for P \Aqt\t=n we need to calculate (c^ ■ ■ ■ Cij k ) c . To do so 
we introduce the generating function 

Z N (h) = Y,W N (C)l[e h ^ , (7) 

C i<j 

with generating fields h = {hij}. This allows us to write 

1 d k 



\ c abi ' ' ' c ab k )c rv 



Z N dh nh , ■ ■ ■ dh„ 



■Z N (h) 



(8) 



h=0 



with Zn = Zjv(O). To carry out the calculation of Z^^h) we introduce a Fourier 
representation for the Kronecker delta constraints appearing in the definition of the 
weight (Ej): 



71 dw-i 

r-r exp I IWa ■ 



L N \ 

t=l j=l J 



(9) 



with u>i = (w^, ...,w^). The sum over C may now be performed explicitly, obtaining 
for large N 



[Cab! 



N . 

Q( k i,kj) ( e 

i, 7=1 ^ 



' TV 

n 



diVi 



x exp 



c 
2iV 



L(2tt) 



x exp 



L 



w 



bi 



(10) 

t=l t=l / _ 

To integrate out the ^-variables and apply saddle-point integration we first need to 
decouple terms comprising vertex indices. To achieve this, we introduce the order 
parameter 



1 N / L 

[k, q) = jjYl 5k *> k exp ( ~ l Yl 
i=i \ i=i 



(11) 
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allowing us to write the average (c^ • • • c a t, k ) c as follows 



{Cabl ' "° abk)c X Y„ qka Hka, <Zi) ■ • -tKfca, q ka )h a ({q^U) 

k * J2 qi ,..., qkh _ l ip{kb j ,qi)---ip{k bj ,qk b -i)Ik b . [K,{Qj}j . , , 
xJI ^ 1 — L\ , (12) 

i=i i } ( k b j ,qi)---^(h j ,qk bj )h bj [{qj}j=i) I 

where we have taken A; = k a , and have introduced the indicator function 

L 

n*({ft}w=n ,5 *w,E*-x^- i) ' (i3) 

which enforces the relationship between the generalised degrees of the neighbours of a 
given vertex, as noted in (j2J). The measure in (fT2|) is defined by 

/{#^}e^)(---) 

V "'*(*,*) f {d<j)di>}e N *(<t>>^ ' 1 J 

with 

$(0, = -c ^ q)<P( k , q) + \ Yl Q ( k > q ^ k > q ^ q > fc ) 

k,q k,q 

+ Y,p(k)\n ^qi)--^{k,q k )h({qt}U) ■ ( 15 ) 

fc q!,...,q k 

In the limit N ^ oo this measure converges to a functional Dirac delta centred at the 
saddle point of $(0,-0). Extremising $(0,-0), we find saddle-point equations 

ip(k,q) = Q(q,k)(f)(q,k) , (16a) 

«*, ,) = p W ^ . (16i) 

; ; E„ B •■■*(*,»)n»({*}J.i) 

Finally, returning to the definition of P ({<fe}t=i, using the preceding result together 
with the saddle-point equations (j!6aj) and (11661) . we reach 

fc N 

P({qt}Li,k) = j^dsr^ 6 ^ Yl 5 k n , qi -5 k]k ,q k {c ih ---c i]k ) c 

»=i h<-<ik 

- P(k) ^ ^ " ' ^ fc ' gfc ^ fc (i^l^lj (17) 

E qil ..., qfe ^(fe, gi) • • • i>{k, q k )ik {{qt}t=i) ' 

It is not our intention to quantify ( |T7|) as this requires the solution of the saddle-point 
equations, which is a cumbersome task for general values L. It is nonetheless interesting 
to keep in mind this result as it will help us to understand the later results for the spectral 
density. The following marginals will also be relevant for our subsequent discussion: 
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(i) The probability P(q, k) of finding a pair of connected vertices with generalised 
degrees k and q: 

P(q, k) = lirn^ (-^Yl 6k ^' k c v 

= Yl p ({vt}t.=i,q^)=Q(k,qMk,qMq,k), (is) 

91 , — ,Qk-l 

and its conditional distribution 

P(q\k) = P(q, k)/P{k) = Q(k, q)<t>{k, q)<f>(q, k) (~p(k)) . (19) 

(ii) The conditional distribution P ({q i }^ =1 |fc): 
P {{q t } k t=1 \k) =P{{q t }t 1 ,k)/P(k) 

i/j{k, q x )--- i/>(k, q k )l k {{q t }t=i) 



Hk, qx) ■ ■ -^(k, q k )I k ({qtjti) ' 
(iii) The conditional distribution P ({qt}t=i\Qi k): 
P {{q t }t!\q, k)=P ({q t } k t=1 , k) /P(q, k) 

tj)(k, q x ) ■ ■ ■ ip(k, q k -x)U (q, {<Zt}t=i) 



(20) 



D*,...,^ ^ 9i) • • • *l>{k, {q, {qt}tJ) 

3.2. Case L = 1 



(21) 



It is instructive to consider the particular case of L = 1 as it has also been discussed 
in previous works (see, for instance, [21, 23J). In this case the generalised degree k is 
simply the degree k. As we only constrain the degrees, the indicator functions in (116&I) 
are identically one, giving 

0(fc,g)= P(fc) il>(k,q) = Q(k,q)<f>(q). (22) 

with P(k) = p(k)k/c being the degree distribution of the nearest neighbour of a vertex. 
Since the right hand side of the first preceding equation is independent of q, so is (f>(k, q), 
and therefore we write 4>(k, q) = (fi(k). Besides, introducing ip{k) = J2 g <?), we reach 

m = Wj> and w = E^'^)- ( 23 ) 

From this, we obtain simplified expressions for ([18]) and ( [191) : 
Moreover, this also induces factorisation in ( 1201 and ( 12TT) . giving 

k 

P{{q t }t 1 \k)=l[P(q t \k), and P ({q t } k ^\q, k) = P {{q t }^\k) . (25) 
t=x 
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Suppose further that we have a separable function Q(k, q) = v(k)v(q). Then, from (J23D 
is easy to see that J2 q v ( ( l) ( P( < l) — ±1j giving the solution ip(k) = ±v(k). This implies 
in turn that P(q\k) = P(q) as expected. 

Having gained some understanding of the typical order parameters involved in the 
problem, we move on to the explicit calculation of the spectral density. 



4. The spectral density of constrained graphs 



Suppose C is the connectivity matrix of a random graph belonging to the ensemble 
under consideration. It is real and symmetric and hence has N real eigenvalues, which 
we denote by {Xf}^ =1 . The natural object of study is the mean spectral density in the 
limit N — > oo, 



P(A) 



/ 1 N 
lim -V5(A 

\ i=l 




(26) 



To compute ( 1261) . we first recast the problem in terms of a Gaussian integral. Following 
Edwards and Jones [27J, one may write 



p(X) = — lim lim 







, r Irn -— (\nZc(\ £ )) r , 



where \ e 



A — ie and Z c is given by 

N 



z c {K 



\\dxi 



i=l 



exp 




E 



C. . rf . fy ■ 



(27) 



(2f 



Ignoring the imaginary units, this object is reminiscent of the partition function 
of a system of dynamical variables interacting on a graph, much like the Gaussian 
ferromagnetic model introduced in [29J. 

Reasoning along these lines, the tools of statistical mechanics can be brought to bear on 
the calculation of spectral density. In particular, the replica method has been frequently 
applied, leading either to approximative schemes [7J, QUI dU [12], or more recent exact 
solutions [El [151 126]. We proceed with the analysis for our ensemble, paying close 
attention to the impact of constraining the generalised degrees. 



4-1- General calculation 

To evaluate the average of the logarithm in (J27j) . we apply the replica method, writing 



(lnZ c (A £ )) c = lim-ln(^(A £ )) c 
where the replicated partition function reads 



c 



N 



Y\_ dx i 



i=l 



Ae 



iY 



ex p -^J2 X ^ 



i=l 




exp 



E 

Kj 



(29) 



(30) 



c 
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with Xi 



x 



(i) 



% i • • • i 



X 



and dxi = dx\ ■ ■ ■dx\ a> . As in the previous calculation, we 



» 



use a Fourier representation for the Kronecker deltas in the weight ([3]) to enable us to 
perform the ensemble average explicitly, writing, for large N, 



{z n c {K)) 



x exp 



c 



1 

Zn 

N 



N 



i=i 



' N 

n 



dw; 



L(2vr) 



1,5=1 ^ ' 



To decouple vertices we introduce the following order parameter 



1 N 

<$>(x,k,q) = — ^25 ki)k S( 

i=i 



[x, - x)c:\}> j — i wf 1 - ) 



(31) 



(32) 



with x = (x^\ . . . , x^) and 5(xi — x) = 5(x^ — x^) ■ ■ ■ 5(x^ — x^). After enforcing 
the order parameter using a functional Dirac delta and rearranging terms we can write 
a compact expression for the replicated partition function: 



(33) 



where 



F($,q) = -cJ2 J dx^{x > k,q)^{x > k,q)-^^2 P (k)p(q)Q(k,q) 

k,q k,q 

+ |E ^ / dxd V^ <?)$(?/, 9- k Y XV (34) 

fc,q 

c k r a 2 

+ ^P(*0 In - / ctee-^* *(a;, fc, 9l ) • • • fc, q k )l k {{q t } k t= i) ■ 



Ql,--,Qk 

The integral (1331) can now be evaluated by steepest descent. Extremising T with respect 
to $ and \1/ we obtain the saddle-point equations 



*(x,k,q) = Q(k,q) dy$(y,q,k)e ix y 



(35 a) 



$(a;,fc,qr) = P(fe)— — > _ — — ; ; — ; - - - • (356) 



fdye~^E c 



^ fc ' 9i) ■ ■ ■ *(Vi 9* ({Qt}ti) 
The natural next step in the calculation is to make a replica symmetric ansatz. In the 
case of graphs with unconstrained topologies, the correct form for the order parameters 
has been recently established as a superposition of Gaussians [T31 [HI [15] . With a modest 
amount of foresight we extend this to the correlated case by writing 



_JL(a;(a0y2 
g 2A \ ' 



$(*, fe, q) = <j>(k, q) J dA 0( A|fc, <?) J] — -==- 

a=l 



*(aj,fe,g)=V(fe,g) / dA^(A|fe,g)J] 



" e -#(^ (a) ) 2 



(36a) 
(366) 
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where we assume that the densities ip(A\k,q) and <f>(A\k,q) are normalised, i.e. 
J dAip(A\k, q) = 1 and similarly for <f)(A\k, q). Note that the parameter A is generally 
a complex variable and dA = aTieAdlmA. Plugging the ansatze into the saddle-point 
equations and taking the replica limit n->0we obtain 

if>(A\k, q)*P(k, q) = Q(k, q)(j>{A\q, k)<j>{q, k) , (37a) 

P(k) 1>{k, q x ) ■ ■ ■ ip(k, qk-i)I h (q, {q t }t=i) 



4>(A\k,q) 



E 



qi, 



,Qk-l 



if>(k, mi) • • ■ if>(k, m fe )I fc ({mt}y 



jfc-1 



Y[dA t ^(A t \k,q t 



t=i 



5 A 



1 



Integrating with respect to A in ( 13 7 d) and (13751) reveals 
^(fc,q) = Q(k,q)(j)(q,k) , 

0(fe,g) = P(fc)- 



(375) 

(38 a) 
(386) 



E mi ,... )m , ^(fc, mi) • - • m fc )Ifc ({m t } t y 

and we recognise these equations as precisely the saddle-point equations (116 al) and (11651) 
in the earlier calculation of the asymptotic graph statistics. This fact is not surprising 
since these objects depend only upon the generalised degrees (k, q), and in this regard 
the order parameters in both calculations are the same. In turn, this also reveals 
that the expression between braces in (1376)) is precisely the conditional distribution 
P ({q t }t=i\q,k). All in all, we can write a self-consistency equation for the density 
Tp(A\q,k): 

^(A\q,k)= P({Qt} k t =?\q,k) 



'k-l 



Y[dA t t{j(A t \k,q t 



t=l 



5 A- 



+ T,U A t 



(39) 



The simple form of this equation is possible thanks to the link between the order 
parameters in this calculation and the asymptotic form of the conditional distribution 
(|2ip . Moreover, written in this way, the physical interpretation is clear: ip(A\q,k) is 
the conditional density of the parameter A for a vertex of generalised degree k, given 
that it has a neighbour of generalised degree q. 

All that remains now is to compute the spectral density. From (127)) we can write 



P(A) 



2 d 
- lim lim — Im — JF($, vp) 

lim -ReVp(fe) / dA ^ phys (A|fc) A 

e^0+ 7T ^-^ J 



with 



ip phys (A\k) 



- E 

91, ■■■!<?& 



P({Qt}Uk) 



Y[dA t i)(A t \k,q t 



t=i 



5 A 



iX e + Et=l A t 



(40) 



• (41) 
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The calculation is now complete. We have expressed the limiting mean spectral density 
for the connectivity matrices of graphs from the ensemble defined by ([3]) in terms of 
densities ^(A|qr, fc), moreover, we have found that the consistency equation ( |39i) for 
these densities is phrased simply in terms of the asymptotic generalised degree statistics 
of the graphs. 

Although, as this calculation has shown, the spectral density of the ensemble studied 
here is governed entirely by its generalised degree statistics, we should point out that 
this is by no means the general rule. 



4.2. Case L = 1 



Let us again consider the particular case L — 1. As before the generalised degree k is 
simply the degree k, and we have P ({q t }^~l\q, k) = P ({qt}t=i\k>) ■ This implies that 
the left hand side of ( 1391) does not depend on q, that is i/;(A\q, k) = i/;(A\k). This 
implies in turn on the right hand side of (1391) that ip(A t \k,q t ) = ip(A t \q t ). Moreover, 
since P ({qt}tZi\k) = Ylt=i Piftlk) anc ^ u P on defining (p(A\k) = J2 q P(q\k)4>(A\q), we 
can write the following self-consistency equation for <p(A\k): 

g-i 

Y[dA t <p(A t \q) <) ( A 



v (A\k) = J2P(i\k) I 

q J 
The spectral density is given by 

P(A) 



(42) 



lim — Re^^p(fc) 



+0+ 7T 



dA V>ph ys (A|/c) A 



with 



^ phys (A|A;) 



l[dA t( p(A t \k) 



t=i 



5 A- 



(43) 



(44) 



The situation is further simplified if Q is taken to be separable, in which case P(q\k) 
P{q) and we recover the results of [H] for the uncorrelated case. 



4-3. Graphs with community structure 



We turn our attention to graphs with communities whose ensemble weight is given by 
(j4j). In this situation it is again possible to compute expressions for the mean spectral 
density of the ensemble using the replica method. Let us briefly outline the main features 
of the calculation. 

Introducing N vectors {Xi)f =ll of M components each, X-i = (x^i, . . . , x^m), one may 
write 



where 

Z C {\e) 



2 d I 

p(X) = — lim lim — — — Im — ( lnZ c (A e ) 

■)c denotes the ensemble average and 



N 



HdXi 



i=l 



exp 



N 



(45) 



(46) 



i=i 



i<j 
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and dXi = Y\ M =1 dxi^ m . We consider the thermodynamic limit to be given by iV — > oo, 
whilst M remains fixed and finite. Treating as individual vector-valued dynamical 
variables the calculation proceeds as usual via the replica method. The order parameters 
take the same form as in the previous calculation, however with the replica symmetric 
ansatz parameterised by M x M matrices A. In the end, one obtains the following 
self-consistency equation: 

"fc-i 

Y[dA t i>(A t \k,q t ,)dB t ii(B t 



i/>{A\q,k) = P({Qt} k t=l\q,k) J 

j dAu{A)5 (a - (^i(\ £ I M - A) + B tA t B, 

The spectral density is recovered via 

p(A) = lim VVfc) / dA% Yiys {A\k)TiA 

e->0+ 7lM * — ' J 
k 

where 



Ql-Qk 



Y[dA t ij(A t \k,q t )dB tf i(B t ) 



J dAu(A)5 f A - U\ £ I M - A )+J2 B t A t B. 



(47) 



(4f 



(49) 



Note that in the case M — 1 (that is, when each community is a single vertex) we have 
v{A) = Sa,o, and /i(-B) = 8b,\ and we recover the result of the previous calculation. Also, 
for general M, taking L = 1 induces the same simplifications as observed previously. 



5. Numerics 



In general the self-consistency equation (j3HD is not exactly solvable. However a numerical 
solution may be efficiently obtained using population dynamics [30| [14"] . Although we 
will consider mainly the case L — 1, we describe this numerical procedure for the general 
case. First, for all possible pairs of generalised degrees (q,k), each density ip(A\q,k) 
in fl39l) is represented by a population of Af variables {Aj(qr, k)}^L v The following 
procedure is then repeated a predefined number of iteration steps: 

1. Choose a pair of generalised degrees (q,k) and a variable A a (q,k) uniformly at 
random from its population. 

2. Randomly select set of generalised degrees {q}\Zi according to the distribution 
P{{q t )U\Q,k). 

3. Choose k — 1 variables {A^ 1 (fe, q\), . . . , A£ fc _ 1 (fe, qk-i)} uniformly at random from 
their populations. 

4. Assign {t\ £ + £ t fi A*(fc, qt))- 1 -> A a (q, k) 
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This procedure is adapted straightforwardly to find a numerical solution of -?/>phys(A|fc) 
using (jHJ, which is then used to calculate the spectral density from fj40|) . 
To assess our results we compare our analytical findings with results from numerical 
diagonalisation of graphs. Although the starting point of our calculations was the 
ensemble definition ([3]), the final equations are phrased only in terms of the resulting 
generalised degree statistics. It is therefore appropriate to compare our results directly to 
data coming from random graphs with prescribed generalised degree statistics, without 
concerning ourselves with intermediate step of determining a suitable choice of weight 
Wn(C). For L = 1 we heuristically adapt the Steger and Wormald algorithm [31] 
[JJ to generate graphs with a given a connected degree-degree distribution P(k, k') in 
the following way: given a degree sequence k = (ki, . . . , k^) with number of edges 
m = | X/i=i iterate the following procedure: 

1. Let E be a set of assigned edges, k = (k\, . . . , k^) an iV-tuple of integers. 

2. Initialise E — 0,k — k 

3. Choose two vertices Vi,Vj G V with probability 

Pij oc P(ki^ kj^jk{kj and (vi,Vj) $ E. 

Reduce ki, kj by 1. 

4. Repeat Step 3 until no more edges can be added to E. 

5. If \E\ < m report failure otherwise output graph. 

The input for this algorithm is the degree sequence of the graphs to be generated, 
however, our results are expressed in terms of the degree distribution p(k). We therefore 
need to generate degree sequences which are compatible with p(k). We discuss two 
possibilities: 

Random degree sequence: For each instance, the degrees are randomly drawn from 
p(k). There is a chance that no graph can be generated exactly fitting the resulting 
degree sequence k. In this case the degree sequence is said to be non-graphical. To 
deal with this, one may choose either to check the graphicality of k before generating 
graphs, or simply accept all graphs generated regardless of whether they fail step 
5. There are various ways to check graphicality, for instance, a theorem of Erdos 
and Gallai states that a degree sequence k with k\ > ■ ■ ■ > k N and J2f=i ^ even is 
graphical if and only if, for all n = 1, iV — 1 

n N 

ki < n(n — 1) + min{fcj, n} . (50) 

i=l i=n+l 

For a discussion on graphicality and the generation of random graphs, see [33J. 

Fixed degree sequence: Select a set of positive integers {N, N±, N2, ■ ■ ■} such that 
iV = Ylk-^k and p(k) ~ N^/N. We then generate random degree sequences with 
Ni vertices of degree one, N 2 vertices of degree 2, and so on. 

I We have also tried an adapted version of the algorithm suggested in [32] with similar results. 
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In our experience, the adapted Steger-Wormald algorithm yields graphs with the desired 
properties and produces almost no failures, provided one selects the appropriate method 
of generating degree sequences. 

5.1. Case L — 1. Correlated degrees 

Let us consider a graph ensemble in which the degrees of neighbouring vertices are 
correlated, that is, P(k, k') does not factorise. For the sake of simplicity, we consider 
graphs whose vertices can only have degrees 2, 3 and 4, with the following degree 
distribution p(k) and conditional distribution P(k'\k): 

18, 4 r 15, 

P[k) = — d k ,2 + ^4,3 + > 

2 4 11 

P{k'\k) = -5 k ,2<5k',2 + <5fc,3#fc',3 + -5k,4$k',4 + -T<5fc,4<5fc',2 + ~8 k ,24',4 • (51) 

3 5 5 3 

To compute the spectral density of the resulting ensemble, we solve the self-consistency 
equation (j4"2l using population dynamics as described previously In this particular 
case, for each value of k G {2,3,4} the corresponding density <p(A\k) is represented by 
a population of TV = 10 4 variables {Ai(k)}fL±, which are iterated over 200 MC steps. 
The spectral density is then computed via (T441) and (pf3l) . To obtain smooth results, the 
spectral density is averaged over a further 50 MC steps. 

For comparison, we have also calculated the spectral density by numerically 
diagonalising 1000 graphs of size N = 2000. In this case, each instance is produced by 
first generating a degree sequence k = (ki, . . . , k N ) according to p(k) and then applying 
the adapted Steger-Wormald algorithm as described previously. We have checked that 
both methods of generating the degree sequences produce equivalent results for large 
samples and that the number of failures of the algorithm is negligible compared to the 
sample size. 

The results of population dynamics and numerical diagonalisation are presented in 
Figure [1] which, apart from peaks at A = 3 and A ~ 3.7 due to finite size effects, 
shows excellent agreement. 

Note from the choice of P(k'\k) that there is a bias towards edges between vertices of the 
same degree, which suggests that the spectral density should share some features with 
the spectral densities of regular graphs of degrees 2, 3 and 4. Indeed, there are peaks 
close to ±2 and ±v^8, coming from peaks in the spectral density of regular graphs of 
degree 2 and 3, and the domain of the spectral density is approximately [ — 2^/3, 2\/3], 
the domain of the spectral density of the regular graph of degree 4. 
We consider next an ensemble with a power-law degree distribution and correlated 
degrees. Such models arise often in the study of real world complex networks. We have 
chosen P(k,k') oc r kk ' , where r < 1, with a maximum value for k being k mSLX . In the 
limits fc max — > oc and r — > 1, this choice results in a power-law degree distribution with 
exponent 2, though for the purpose of simulations, we will take r = 0.999, and keep 
^max finite. 

We use the previously explained algorithm to generate graphs with distribution P(k, k'). 
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Figure 1: Comparison of the results of population dynamics (blue line) and direct 
diagonalisation (red histogram) for the choice of P(k'\k) given in (|5"T|) . To construct the 
histogram we use the adapted Steger-Wormald algorithm to generate 1000 graphs of 
size N = 2000. The degree sequences were generated randomly according to p(k) given 
in (l5Tj) and only eight failures were reported. 

For graphs of size N = 2000 it is necessary to take a rather low maximum degree 
fc max = 45 ( ~ VN if tax is taken any larger, additional correlations occur between 
high degree vertices [34], and the failure rate becomes unacceptable. Alternatively, we 
could have larger values of k max by increasing the graph size, but that would make 
the numerical diagonalisation computationally expensive. The degree sequences are 
generated randomly and checked for graphicality before generating the graph. 
The spectral density can again be computed via population dynamics. Figure [2] shows a 
comparison between the results of population dynamics and a histogram of eigenvalues 
from 500 random graphs of size N = 2000. In this case we take a population of M = 10 3 
variables {A^/c)}^, which are iterated over 100 MC steps. The spectral density is then 
averaged over a further 50 MC steps. 

As before some of the salient features of the spectral density can be intuitively explained 
in terms of the underlying graph structure. First we note that the spectral density for 
fcmax = 45 and r = 0.999 has a bounded support as expected due to the Perron-Frobenius 
theorem, but as we take the /c max — > oo and r — > 1, the mean degree diverges and the 
spectral density presents heavy tails. It is interesting to analyse the contribution of 
vertices of high degree to the spectral density. As mentioned in [8], for such vertices 
it may be sufficient to consider only the mean behaviour of the neighbouring vertices. 
In this effective medium approximation (EMA) the approximate behaviour of the tails 
for very large |A| reads p(A) ~ 2k\p(k\)/\\\, where k\ = A 2 + 0(1). A more rigorous 
analysis in [35J states that the largest eigenvalues of graphs with heavy-tailed degree 
distributions occur close to the square roots of the largest degrees. 
In Figure [3] we show the results of population dynamics simulations in the tail of the 
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Figure 2: Comparison of the results of population dynamics (blue line) and direct 
diagonalisation (red histogram) for the choice P(k,k') oc r kk , where r = 0.999 and 
fc max = 45. The adapted Steger-Wormald algorithm was used to generate 500 graphs of 
size N = 2000, whose eigenvalues were used to construct the histogram. 

spectral density for a much larger maximum degree of fc max = 400. The approximate 
curve given by the EMA gives a reasonable fit with the result of the simulation and, as 
expected, the density drops dramatically shortly after y/k max = 20. The contributions 
to the density coming from high degree vertices can be isolated in the output of the 
population dynamics algorithm; we have included in Figure [3] contributions from several 
high degrees k, each of which exhibits sharp peak close to \fk. 

The other main feature of the density shown in Figure [2] is the presence of Dirac delta 
peaks at -1, and 1 whose weight may be bounded by using the distributions p(k) and 
P(k\k'). For instance, the weight to the peaks at ±1 has contributions from connected 
pair of vertices of degree 1, which for the choice fc max = 45 and r = 0.999 have a 
likelihood of p(l)P(l|l)/2 = 0.0031, not far from the exact weight of 0.0037 obtained 
from numerical diagonalisation. A similar intuitive argument can be used to obtain 
a bound for the weight of the Dirac delta peak at zero, whose appearance is due to 
dead-end vertices [H [36] . 

5.2. Case L = 2. Levels of approximation 

Suppose that the exact knowledge of a graph ensemble is reduced solely to a set of 
statistical properties captured by, for instance, the degree distribution p{k) and the 
conditional distribution P{k\k'). Whilst in a few cases such quantities suffice to fully 
characterise the graph ensemble, this is not generally true. We would like to understand 
in which way the lack of more accurate information affects the spectral density. 
With this in mind let us consider a graph ensemble with generalised degrees {ki\f =1 
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Figure 3: Detail of the tail of the spectral density for the choice P(k, k 1 ) oc r kk> , where 
r = 0.999 and fc max = 400. The continuous blue line is the full result from population 
dynamics, with the labelled black lines being isolated contributions from vertices of high 
degree. The dashed red line shows an estimate for the tail derived from the effective 
medium approximation. 



such that ki G {k 2 , ^3, ^4}, where 

«3 = [ I ) , «a = I ; I , ; I I . (52) 





This ensemble is composed of graphs with vertices of degrees 2, 3 and 4. Moreover, 
those vertices of degree 2 must be connected to one vertex of degree 3 and one of degree 
4 and those of degrees 3 and 4 must be connected to vertices of degree 3 only. A portion 
of such graph is shown in the leftmost part of Figure HI 

For this graph ensemble a quick counting argument gives the following expressions for 
p(k) and P(k\k') 

12 4 3 

P(k) = — S k ,2 + —^,3 + —4,4 , (53a) 

P(k\k') = S k>2 6k',4 + 4,24',3 + 0^,3 + J , (536) 

knowledge of which does not fully characterise the graph ensemble. Note that this 
example is such that the set of self-consistency equations fl39|) can be solved exactly. To 
do so we first need the conditional distribution P ({qt}tZi\Q> k) which in this case reads 

P(k 3 |k 4 ,k 2 ) = 1, P{k 4 \k 3 ,k 2 ) = 1 , 

P(k 2 , k,2\k 2 , K3) = 1 , P(k 2 , «2, k 2 |k 2 , « 4 ) = 1 , 

or zero otherwise. This results in a set of self-consistency equations for the densities 
{■0(A|«3, k 2 ), i/j(A\k 2 , k 3 ), ^(A|/«4, k 2 ), i/;(A\k 2 , k^)} that admits a solution of the type 
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Figure 4: Left - a typical neighbourhood of a vertex of degree 3 of a random graph 
specified by (1521) . Middle - a neighbourhood of a vertex of degree 3 in a random graph 
with degree distribution (153 aft and P(k\k') given by (15361) . Right - a neighbourhood of 
a vertex of degree 3 in an uncorrelated random graph with degree distribution ( 153 al) . 



i[)(A\K a , K b ) = 5(A — A a b ) with A afe obeying a simple set of algebraic equations 
1 1 1 1 



A 



3.2 



A 



v 2,3 



»2,4 



zA £ + A 2 , 4 *A e + A 2i3 i\ £ + 2A 3i2 «A e + 3A 4j2 

Finally, to find an expression for the spectral density we first need the distribution 
f | fe) which in this case reads 

1 1 



., ■ ^(«4,«3|«2) ,, 



P(« 3 , K 4 |«2) 
P(« 2 , «2, «2|«3) = 1 , P(«2, «2, «2, «2 1 «4) = 1 

or zero otherwise. This yields 



p(A) = lim Re 

Hy ' £ ^Q+ 197T 



12 



I | (55) 

_i\ £ + A 2 , 3 + A 2 , 4 «A £ + 3A 3i2 iA £ + 4A 4 , 2 J 

Upon solving (I54l) . plugging the solutions into (|55|) and carefully analysing the poles, 
we can write 



P(A) 



^(A) + U(\ + V3) + ±6(A - VS) 



12|2A 2 - 7| V-25 - A 2 (-7 + A 2 ) (14 - 7A 2 + A 4 ) _ , , 

H ~ i wxo TTTTo ^TTTn ^Ti %(Aj 



19tt|A(A 2 -4)(A 2 -7)(A 2 -3)| 
with ft = [A;, Al] U [At, \X\ and 



A^ 



14 + 2/iV21 + 8a\/6 



(56) 



(57) 



with a, \i G {— , +} and where I^(A) = 1 if |A| G 1Z or zero otherwise. 

This example was chosen specifically to keep the local structure of the graphs 

deterministic and thus make equations (1391) exactly solvable in the case L = 2. If 
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we reduce the statistical information we have about the ensemble, either by taking 
L — 1, or assuming the degrees to be uncorrelated, the local graph structure becomes 
random, and the resulting expressions are no longer exactly solvable. Figure H] shows 
the randomising effect of these simplifications. 

We have used population dynamics to compute the spectral density of random 
graphs with degree distribution (153 and degree-degree correlations (153 fej) . as well 
as uncorrelated random graphs with degree distribution (153 a\i . For each degree, the 
population size is TV = 10 4 , iterated over 200 MC steps. The spectral density is then 
averaged over 50 MC steps. The results are shown in Figures Mjo) and[5](c), alongside 
histograms of the eigenvalues of 1000 graphs of size N = 1900, generated using the 
adapted Steger-Wormald algorithm 0. Figure 0(a) shows the exact spectral density 
given by ( l56l) . plotted alongside a histogram of eigenvalues obtained from randomly 
generated graphs of that type. 

In Figure [5], the effects of reducing knowledge of an ensemble (and hence increasing 
randomness) are clearly visible. In addition to a general smoothing effect, which one 
might expect, the most striking feature is the appearance of gaps in the spectral density. 
When the ensemble is fully specified by ( 1521) . the continuous part of the spectral density 
is divided into four disjoint components. When one specifies only the degree distribution 
and degree-degree correlations, the number of components reduces to two, and when only 
the degree distribution is known, there is no gap in the density at all. 
The appearance of the gap can be traced back to the periodicity in the random graph 
ensembles (see [37]). For instance, in the original ensemble a walker moving away from a 
central vertex will repeatedly visit vertices of degree sequence {...,2,4, 2, 3, 2, 4, 2, 3, . . .} 
of periodicity four which is also manifestly explicit in the equations (J54l) . In the case that 
the degree-degree correlations are specified, whilst the degree sequence of a walk is now 
random and therefore not strictly periodic, it is still true that every other vertex visited 
will have degree 2 (see Figure HI), periodicity which is enough to split the spectral density 
into two components. In the last case, where only the degree distribution is known, the 
sequence is fully random (see Figure HI), and the resulting spectral density has no gap. 



5.3. Community structure 

The population dynamics algorithm used to solve (l39l is easily adapted to solve 
the equivalent self-consistency equation ( HTI) for the communities model; one simply 
initialises populations of M x M matrices A. and updates them according to (I4"T|) . 
The presence of communities in a graph typically results in a very different spectral 
density, to illustrate this we consider a simple choice for the community structure 
ensemble. Suppose we have communities given by the complete graph on M vertices, 
connected in a Poissonian random graph of average degree c, in which connected 
communities are joined by a single randomly drawn edge. In the weight (jl]), this 

% It is necessary to slightly modify the adapted Steger-Wormald algorithm to make sure the constraints 
specified by ([52]) are met. 
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(a) Comparison of the exact spectral density given by equation ([55)1 (blue line) and direct 
diagonalisation (red histogram). We have checked the isolated peaks at A ~ ±v7 are 
due to finite size effects. 
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(b) Comparison of population dynamics with degree distribution (|53q|) and degree-degree 
correlations (|536[) (blue line) and direct diagonalisation (red histogram). To visualise the 
Dirac delta peak we have taken e — 10~ 6 . 
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(c) Comparison of population dynamics with degree distribution only (blue line) and 
direct diagonalisation (red histogram). 

Figure 5: Comparison theoretical results (blue lines) and direct diagonalisation (red 
histograms) for the different levels of approximation to the ensemble specified by ( 152]) . 
In plot (a), the blue line shows the exact result for the spectral density, given by (156|) . 
in plots (b) and (c) the result of population dynamics is shown. We use the method 
outlined in the text to generate graphs for plot (a), and the adapted Steger-Wormald 
algorithm for (b) and (c). In each case, 1000 graphs of size iV = 1900 were generated 
and diagonalised. 
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Figure 6: Comparison of the results of population dynamics (blue line) and direct 
diagonalisation (red histogram) for the community structure ensemble described in the 
text. The grey curve shows the high connectivity limit for this ensemble. 



1 for all v(A) 



na^a,o)(n a ^&^w)> 



corresponds to the choices Q(ki,kj) 
and fi(B) = M- 2 J2 a , b $B ab ,i U( c ,d)^(a,b) S B cd ,o ■ 

Taking M = 5 and c = 5, we use population dynamics to solve (1471) for this ensemble, 
the spectral density is them computed using fj4"8]) . To compare with the results of direct 
diagonalisation, 1000 graphs of size N = 5000 were generated. A histogram of their 
eigenvalues, alongside the result of population dynamics, is shown in Figure [6j 
It is well known that in the limit c — > oo the spectral density of a Poissonian 
random graph with average degree c, and edges of weight 1/a/c, converges to Wigner's 
semi-circular distribution [5j [13]. We can compute a generalisation of this result 
for the community ensemble considered here through an appropriate treatment of 
the self-consistency equation (j47|) . It is necessary to re- weight the edges between 
communities, in order to keep the spectral density bounded as c — ► oo, we take 



MS) = M- 2 E a ^ B 



'M/c 



ll(c dWa b) ^B cd ,o ■ Keeping only the terms relevant in the 



oo limit we obtain an expression for the mean A. 
(A) = (i(X e I M - K M ) + c(b(A)B t 



(5f 



where Km is the connectivity matrix of the complete graph on M vertices. For the 
above choice of fx(B), we have (BAB T ) = (l/cM)AI M , where A = (l/M)Tt(A). 
Diagonalising K M , we obtain a cubic equation for A, 



A 



In the case M 



M - 1 

M 

5, we solve 



1 



1 

M 



1 



i\ £ + A + i ' \Mji\ £ + A-(M 
59]) to find the following expression: 
A 2 -3A + 18 



1 i 



u 



9u 



MA) 



(59) 



(60) 
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where 



-4 A 3 - -A 2 - 2A + 4 + — V3d 
27 6 18 



1/3 



u = 



(61) 



and d = -25 A 4 + 156A 3 + 72 A 2 - 1296A + 864. The domain is given by I^(A) = 1 if 
d > 0, and zero otherwise. 

6. Conclusions 

Past work on the spectral density of random graphs has typically been confined to 
simple ensembles in which at most the degree distribution is specified and even then 
exact results have only been obtained relatively recently. At the same time, the field 
of complex networks has gained a great deal of attention from all over the scientific 
community. Before the analysis of spectral density can be used to provide insights into 
the behaviour and characterisation of complex networks, progress must first be made to 
expand the class of ensembles for which exact results are obtainable. 
In this paper we have sought to do just this, through the calculation of the spectral 
density of random graphs with constrained topologies. Complex correlations between 
the degrees of non-neighbouring vertices are incorporated in the constrained generalised 
degree ensemble and we also introduce a simple extension of this model to one featuring 
a community structure. The important statistical properties of the constrained degree 
ensemble are captured in the distribution (jSJ), which we compute via a saddle point 
analysis in the large N limit. This calculation foreshadows the replica calculation of the 
spectral density and provides important insight. 

For the problem of determining the mean spectral density, we take standard steps to 
map the problem onto one of an interacting system of dynamical variables, to which 
the replica method is applied. Following recent advances, the form of replica symmetric 
ansatz is identified as a superposition of Gaussians. Exploiting the insights gained in 
the earlier calculation, we obtain closed expressions for the spectral density in terms of 
the statistical properties of the graph ensemble fl39|) . Similar equations are found for 
the community structure ensemble ( 1471) . 

Though the resulting equations may not often have easily found analytic solutions, 
they can be efficiently solved numerically using the population dynamics algorithm we 
describe earlier. In this way it is possible to analyse the spectral density of a given graph 
ensemble without the need to generate and diagonalise large numbers of graphs. An 
instance of this is provided by the discussion of the tails in the example with a power- 
law degree distribution; here the results of population dynamics could not feasibly be 
obtained by diagonalising random matrices (for fc max = 400, we would require graphs of 
around 1.6 x 10 5 vertices). We hope that the methods discussed here will prove useful 
tools in the study of complex networks. 

Although, as we have demonstrated, the statistics of generalised degrees can have a 
significant impact on the spectral density of the graph ensemble, this is certainly not 
the only factor at work. One aspect of the topology of complex networks which may 
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play an important role, but which we have not considered so far, is the statistics of the 
loops in the graph. Unfortunately, knowledge of the generalised degrees of a graph gives 
no information about loops, and hence the constrained generalised degree ensemble is 
not likely to be useful to study of their effect on spectral density. Looking to the future, 
it seems the next major step forward in the analysis of spectral density of random 
graphs will require techniques capable of handling the effects of loops. Several new 
techniques to correct for the presence of loops in related problems have a appeared 
recently [38j ESI HO] , and we hope that similar ideas may also be applied to the study 
of spectral density. 
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